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Anisotropies in the diffuse gamma-ray background measured by the Fermi LAT 
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The contribution of unresolved sources to the diffuse gamma^ray background could induce 
anisotropies in this emission on small angular scales. We analyze the angular power spectrum 
of the diffuse emission measured by the Fermi LAT at Galactic latitudes \b\ > 30° in four energy 
bins spanning 1 to 50 GeV. At multipoles £ > 155, corresponding to angular scales < 2°, angular 
power above the photon noise level is detected at > 99.99% CL in the 1-2 GeV, 2-5 GeV, and 5- 
10 GeV energy bins, and at > 99% CL at 10-50 GeV. Within each energy bin the measured angular 
power takes approximately the same value at all multipoles £ > 155, suggesting that it originates 
from the contribution of one or more unclustered source populations. The amplitude of the angular 
power normalized to the mean intensity in each energy bin is consistent with a constant value at all 
energies, Cp /(/) 2 = 9.05 ± 0.84 X 10“ 6 sr, while the energy dependence of Cp is consistent with the 
anisotropy arising from one or more source populations with power-law photon spectra with spectral 
index r s = 2.40 ± 0.07. We discuss the implications of the measured angular power for gamma-ray 
source populations that may provide a contribution to the diffuse gamma-ray background. 
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I. INTRODUCTION 

The origin of the all-sky diffuse gamma^ray emission 
remains one of the outstanding questions in high-energy 
astrophysics. First detected by OSO-3 [1], the isotropic 
gamma-ray background (IGRB) was subsequently mea- 
sured by SAS-2 [2], EGRET [3, 4], and most recently 
by the Large Area Telescope (LAT) onboard the Fermi 
Gamma-ray Space Telescope (Fermi) [5]. The term 
IGRB is used to refer to the observed diffuse gamma- 
ray emission which appears isotropic on large angular 
scales but may contain anisotropies on small angular 
scales. The IGRB describes the collective emission of 
unresolved members of extragalactic source classes and 
Galactic source classes that contribute to the observed 
emission at high latitudes, and gamma-ray photons re- 
sulting from the interactions of ultra-high energy cosmic 
rays with intergalactic photon fields [6] . 

Confirmed gamma^ray source populations with re- 
solved members are guaranteed to contribute to the 
IGRB at some level via the emission from fainter, unre- 
solved members of those source classes. In the EGRET 
era the possibility that blazars are the dominant con- 
tributor to the IGRB intensity was extensively studied 
(e.g., [7-9]), however the level of the blazar contribution 
remains uncertain, with recent results suggesting differ- 
ent energy-dependent contributions from blazars, which 
amount to as little as ~ 15% or as much as ~ 100% of 
the Fermi-measured IGRB intensity, depending on the 
energy [10-13]. Star-forming galaxies [14] and gamma- 
ray millisecond pulsars [15] may also provide a significant 
contribution to the IGRB at some energies. However, 
substantial uncertainties in the properties of even con- 
firmed source populations present a challenge to estimat- 
ing the amount of emission attributable to each source 
class, and currently the possibilitv that the IGRB in- 
cludes an appreciable contribution from unknown or un- 
confirmed gamma-ray sources, such as dark matter anni- 
hilation or decay (e.g., [16-22]), cannot be excluded. 

The Fermi-measured IGRB energy spectrum is rela- 
tively featureless, following a simple power law to good 
approximation over a large energy range (~ 200 MeV to 
~ 100 GeV) [5]. As a result, identifying the contribu- 
tions from individual components based on spectral in- 
formation alone is difficult. However, in addition to the 
energy spectrum and average intensify , the IGRB con- 
tains angular information in the form of fluctuations on 
small angular scales [23]. The statistical properties of 
these small-scale anisotropies may be used to infer the 
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presence of emission from unresolved source populations. 

If some component of the IGRB emission originates 
from an unresolved source population, rather than from 
a perfectly isotropic, smooth source distribution, jbhe dif- 
fuse emission will contain fluctuations on small angular 
scales due to the varying number density of sources in 
different sky directions. Unlike the Poisson fluctuations 
between pixels in a map of a truly isotropic source dis- 
tribution (which we shall call “photon noise” ) , which are 
due to finite event statistics, the fluctuations from an un- 
resolved source population are inherent in the source dis- 
tribution and will not decrease in amplitude even in the 
limit of infinite statistics. Hence, with sufficient statisr- 
tics, these fluctuations could be detected above those ex- 
pected from the photon noise, and could be used to un- 
derstand the origin of the diffuse emission. 

The angular power spectrum of the emission provides a 
metric for characterizing the intensity fluctuations. For a 
source population modeled with a specific spatial and lu- 
minosity distribution, the angular power spectrum can be 
predicted and compared to the measured angular power 
spectrum; in this way an anisotropy measurement has the 
potential to constrain the properties of source popula- 
tions. Other approaches to using anisotropy information 
in the IGRB have also been considered. For example, the 
1-point probability distribution function (PDF), i.e. the 
distribution of the number of counts per pixel, is an alter- 
native metric to characterize the fluctuations [13, 24, 25]. 
In addition, cross-correlating the gamma-ray sky with 
galaxy catalogs or the cosmic microwave background can 
be used to constrain the origin of the emission [26]. 

In recent years theoretical studies have predicted the 
angular power spectrum of the gamma-ray emission from 
several known and proposed source classes. Established 
astrophysical source populations such as blazars [27-29], 
star- forming galaxies [30], and Galactic millisecond pul- 
sars [31] have been considered as possible contributors 
to the anisotropy of the IGRB. In addition, it has 
been shown that the annihilation or decay of dark mat- 
ter in Galactic subhalos [32-34] and extragalactic struc- 
tures [23, 27, 29, 34-39], may generate an anisotropy sig- 
nal in diffuse gammarray emission. Interestingly, the pre- 
dicted angular power spectra of these gamma-ray source 
classes in the multipole range of £ ~ 100-500 are in most 
cases fairly constant in muitipole (except for dark matter 
annihilation and decay signals, e.g. [23, 27, 38]), although 
the amplitude of the predicted anisotropy varies between 
source classes. This multipole-independent signal arises 
from the Poisson term in the angular power spectrum, 
which describes the anisotropy from an unclustered col- 
lection of point sources. The multipole-independence of 
the predicted angular power spectra therefore indicates 
that the expected degree of intrinsic clustering of these 
gamma-ray source populations has a subdominant effect 
on the angular power spectra in this multipole range. The 
angular power spectra of dark matter annihilation and 
decay signals are predicted to be smooth and relatively 
featureless, with the angular power generally falling off 
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more quickly with multipole than Poisson angular power. 

In this work we present a measurement of the angular 
power spectrum of the high- latitude emission detected by 
the Fermi LAT, using ~ 22 months of data. The data 
were processed with the Fermi Science Tools [56], and 
binned into maps covering several energy ranges. Re- 
gions of the sky heavily contaminated by Galactic diffuse 
emission and known point sources were masked, and then 
angular power spectra were calculated on the masked sky 
for each energy bin using the HEALPix package [57], de- 
scribed in [40]. 

To understand the impact of the instrument response 
on the measured angular power spectrum, several tai- 
lored validation studies were performed for this analy- 
sis. The robustness of the anisotropy analysis pipeline 
was tested using a source model with known anisotropy 
properties that was simulated to include the effects of the 
instrument response and processed with the same analy- 
sis pipeline as the data. The data processing was cross- 
checked to exclude the presence of anisotropies created 
by systematics in the instrument exposure calculation by 
using an event-shuffling technique (as used in [41]) that 
does not rely on the Monte- Carlo-based exposure calcu- 
lation implemented in the Science Tools. In addition, val- 
idation studies were performed to characterize the impact 
of foreground contamination, masking, and inaccuracies 
in the assumed point spread function (PSF). 

We use a set of simulated models of the gamma-ray sky 
as a reference, and compare the angular power spectrum 
measured for the data to that of the models to identify 
any significant differences in anisotropy properties. Fi- 
nally, we compare the predicted anisotropy for several 
confirmed and proposed gamma^ray source populations 
to the measured angular power spectrum of the data. 

The data selection and map-making procedure are de- 
scribed in §11, and the angular power spectrum calcu- 
lation is outlined in §IIL The event-shuffling technique 
is presented in §IV, and the details of the models sim- 
ulated to compare with the data are given in §V. The 
results of the angular power spectrum measurement and 
the validation studies are presented in §VI. The energy 
dependence of the anisotropy is discussed in §VII, and the 
implications of the results for specific source populations 
are examined in §VIII. The conclusions are summarized 
in §LK. 


II. DATA SELECTION AND PROCESSING 

The Fermi LAT is designed to operate primarily as a 
survey instrument, featuring both a wide field of view 
(~2.4 sr) and a large effective area (> 7000 cm 2 for 
normally- incident photons above 1 GeV). The telescope 
is equipped with a 4x4 array of modules, each consisting 
of a precision tracker and calorimeter, covered by an anti- 
coincidence detector that allows for rejection of charged 
particle events. Full details of the instrument, including 
technical information about the detector, onboard and 


ground data processing, and mission-oriented support, 
are given in [42] . 

We selected data taken from the beginning of scientific 
operations in early- August 2008 through early- June 2010, 
encompassing over 56.6 Ms of live time [58]. We selected 
only “diffuse” class [42] events to ensure that the events 
are photons with high probability, and restricted our 
analysis to the energy range 1-50 GeV where the PSF of 
the LAT is small enough to allow for sufficient sensitivity 
to anisotropies at small angular scales. The upper limit 
of 50 GeV was chosen because the small photon statis- 
tics above this energy severely limit the sensitivity of the 
analysis at the high multipoles of interest. The data and 
simulations were analyzed with the LAT analysis soft- 
ware Science Tools version v9rl5p4 using the standard 
P6-V3 LAT instrument response functions (IRFs). De- 
tailed documentation of the Science Tools is given in [59] . 

In order to both promote near uniform sky exposure 
and to limit contamination from gamma rays originat- 
ing in Earth’s atmosphere, the tool gtmktime was used 
to remove data taken during any time period when the 
LAT rocked to an angle exceeding 52° with respect to 
the zenith, and during any time period when the LAT 
was not in survey mode. Beginning in its second year of 
operation (September 2009), Fermi has been operating 
in survey mode with a large rocking angle of 50°, in con- 
trast to the 35° rocking angle used during the first year 
of operation. The rocking- angle cut is used to limit the 
amount of contamination from gamma rays produced in 
cosmic-ray interactions in the upper atmosphere by us- 
ing only data taken when the Earth’s limb was outside 
of the field of view (the Earth’s limb has zenith angle 
~ 113°). However, due to the LAT’s large field of view, 
some Earth- limb gamma rays may be observed even when 
the rocking angle constraint is not exceeded, thus the gt- 
select tool was also used to remove each individual event 
with a zenith angle exceeding 105°. We note that all 
events in the data set were detected while the Fermi 
spacecraft was outside of the South Atlantic Anomaly 
region in which the cosmic-ray fluxes at the altitude of 
Fermi are significantly enhanced. 

In order to balance the need for a large effective area 
with the need for high angular resolution, the LAT uses 
a combination of thin tracker regions near the front of 
the instrument and thicker tracker regions in the back of 
the detector. While the effective area of each region is 
comparable, the width of the PSF for events detected in 
the front trackers is approximately half that of events de- 
tected in the back of the instrument. For a measurement 
of the angular power at high multipoles, it is thus nec- 
essary to differentiate between photons observed in the 
front and back trackers of the Fermi LAT. In this study, 
we processed front- and back-converting events sepa- 
rately, using the gtselect tool to isolate each set of events 
and calculating the exposure maps independently. The 
P6_V3_DIFFUSE:FRONT and P6_V3JDIFFUSE:BACK 
IRFs were used to analyze the corresponding sets of 
events. 
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Taking the selection cuts into account, the integrated 
live time was calculated using gtltcube. We chose a pixel 
size of 0.125°, which produces a HEALPix map with res- 
olution parameter A^e = 512. At this resolution, the 
suppression of angular power from the pixel window func- 
tion is subdominant with respect to the suppression from 
the LAT PSF. We adopted an angular step size cos(0) = 
0.025 in order to finely grid the exposure map for differ- 
ent gamma-ray arrival directions in instrument coordi- 
nates. The exposure was then calculated using gtezpcube 
with the same pixel size, for 42 logarithmic energy bins 
spanning 1.04 GeV - 50.0 GeV. These finely-gridded en- 
ergy bins were then summed to build maps covering four 
larger energy bins, as described in §111 A. Using the GaR- 
DiAn package [43] , both the photon counts and exposure 
maps were converted into HEALPix-format maps with 
N eide = 512. 


m. ANGULAR POWER SPECTRUM 
CALCULATION 

We consider the angular power spectrum Ct of an in- 
tensity map I (i!>) where rp denotes the sky direction. 
The angular power spectrum is given by the coefficients 
Ct = (|a* m | 2 ), with the determined by expanding 
the map in spherical harmonics, 

m = Y. a ^yt m {^y a) 

im 

The intensity angular power spectrum indicates the di- 
mensionful size of intensity fluctuations and can be com- 
pared with predictions for source classes whose collective 
intensity is known or assumed (as in, e.g., [31]). The in- 
tensity angular power spectrum of a single source class is 
not in general independent of energy due to the energy- 
dependence of the mean map intensity (/). 

We can also construct the dimensionless fluctuation 
angular power spectrum by dividing the intensity angu- 
lar power spectrum Ct of a map by the mean sky inten- 
sity (outside of the mask, for a masked sky map) squared 
(I) 2 . The fluctuation angular power spectrum character- 
izes the angular distribution of the emission independent 
of the intensity normalization. Its amplitude for a single 
source class is the same in all energy bins if all members 
of the source class share the same observed energy spec- 
trum. since this results in the angular distribution of the 
collective emission being independent of energy. Energy 
dependence in the fluctuation angular power due to vari- 
ation of the energy spectra between individual members 
of the population is discussed in §VII. 

A. Energy dependence 

We calculate the angular power spectrum of the data 
and simulated models in four energy bins. Using mul- 
tiple energy bins increases sensitivity to source popula- 


tions that contribute significantly to the anisotropy in 
a limited energy range, and may also aid in the inter- 
pretation of a measurement in terms of a detection of 
or constraints on specific source populations [39, 44]. In 
addition, the detection of an energy dependence in the 
fluctuation angular power spectrum of the total emission 
(the anisotropy energy spectrum) may be used to infer 
the presence of multiple contributing source classes [45] . 
In the case that a single source population dominates 
the anisotropy over a given energy range, the energy de- 
pendence of the intensity angular power spectrum can 
indicate the energy spectrum of that contributor. 

Since the LAT’s angular resolution and the photon 
statistics depend strongly on energy, the sensitivity of 
the analysis is also energy-dependent: at low energies 
the LAT’s PSF broadens, resulting in reduced sensitiv- 
ity to small-scale anisotropies, while at high energies the 
measurement uncertainties are dominated by low statis- 
tics. We calculate angular power spectra in the energy 
bins 1.04-1.99 GeV, 1.99-5.00 GeV, 5.00-10.4 GeV, and 
10.4-50.0 GeV. The map for each energy bin for the an- 
gular power spectrum analysis was created by summing 
the corresponding maps produced in finely-gridded en- 
ergy bins, as described in §11. 


B. Angular power spectrum of a masked sky 

The focus of this work is to search for anisotropies 
on small angular scales from unresolved source popula- 
tions, hence the regions of the sky used in this analysis 
were selected to minimize the contribution of the Galac- 
tic diffuse emission from cosmic-ray interactions and the 
emission from known sources. A mask excluding Galactic 
latitudes |6| < 30° and a 2° angular radius around each 
source in the 11-month Fermi LAT catalog (1FGL) [46] 
was applied prior to performing the angular power spec- 
trum calculations in all energy bins. The fraction of the 
sky outside of this mask is f B ^ y = 0.325. The 2° angular 
radius for the source masking approximately corresponds 
to the 95% containment angle for events at normal inci- 
dence at 1 GeV (front/back average for P6-V3 IRFs); 
the containment angle decreases with increasing energy. 
The effect of the mask on the angular power spectra is 
discussed below and in §VI F, and the impact of varia- 
tions in the latitude cut is assessed in §VIE. An all-sky 
intensity map of the data in each energy bin is shown in 
Fig. 1, both with and without applying the default mask. 

The angular power spectra of the masked maps were 
calculated using HEALPix, after first removing the 
monopole and dipole terms. To approximately correct 
for the power suppression due to masking, the raw angu- 
lar power spectra output by HEALPix were divided by 
the fraction of the sky outside the mask, f B ^ y . This cor- 
rection is valid at multipoles greater than ~ 10, where 
the power spectrum of the signal varies much more slowly 
than the window function, as detailed below. 

When a fraction of the sky is masked, the mea- 
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FIG. Is All-sky intensity maps of the data in the four energy bins used in this analysis, in Galactic coordinates; the map 
projection is Mollweide. The data shown are the average of the maps of the front- and back-converting events, and are shown 
unmasked (left panels) and with the default mask applied (right panels ). The mask excludes Galactic latitudes |6| < 30° and 
a 2° angular radius around each source in the 1FGL catalog. The map images shown have been downgraded in resolution to 
AT 8ide = 128 to improve the visual quality of the images; however, the analysis was performed on the higher resolution maps as 
described in the text. 
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sured spherical harmonics coefficients are related to the 
true, underlying spherical harmonics coefficients, a£™ e , 
via a matrix multiplication, a* m - 
where is the so-called coupling matrix given by 

Wutfnmt = d 2 n^(n)l^/ m /( n), where the integral 

is done only on the unmasked sky whose solid angle is 
n o b s . Then, HEALPix returns a raw angular power spec- 
trum, C] aw = (2t + l) -1 Ylm Km| 2 , whose ensemble av- 
erage is related to the true power spectrum, Cj rue , as 

{cd = 2ZTI E ct r E 1 w «'— ' i 2 - ( 2 ) 

V mm' 

Now, for a given mask, one may calculate 

Emm and estimate Cf ue by inverting 

this equation. This approach is called the MASTER 
algorithm [47], and has been shown to yield unbiased 
estimates of C| rue . However, while this is an unbiased 
estimator, it is not necessarily a minimum- variance 
one. In particular, when the coupling matrix is nearly 
singular because of, e.g., an excessive amount of mask or 
a complex morphology of mask, this estimator amplifies 
noise. We observed this amplification of noise when 
applying the MASTER algorithm to our data set. 
Therefore, we decided to use an approximate, but less 
noisy alternative. It is easy to show (see, e.g., Eq. A3 
of [4Si) that, when £ TOm , |WV mm '| 2 peaks sharply at 
t = t* and C| rue varies much more slowly than the width 
of this peak, the above equation can be approximated as 

{CD * = cDfsky (3) 

This approximation eliminates the need for a matrix in- 
version. We have verified that this method yields an 
unbiased result with substantially smaller noise than the 
MASTER algorithm at £ > 10. We adopt this method 
throughout this paper. 

C. Window functions 

The angular power spectrum calculated from a map 
is affected by the PSF of the instrument and the pix- 
elization of the map, encoded in the beam window func- 
tion R 7beam and the pixel window function WJ* 1 * respec- 
tively both of which can lead to a multipole-dependent 
suppression of angular power that becomes stronger at 
larger multipoles. Depending upon whether the power 
spectrum originates from signal or noise, corrections for 
the beam and pixel window functions must be applied 
to the measurement differently. For our application, we 
must not apply any corrections to the photon shot noise 
(Poisson noise) term, while we must apply both the beam 
and pixel window function corrections to the signal term 
from, e.g., unresolved sources. While it is obvious why 
one must not apply the beam window correction to the 
photon noise term, it may not be so obvious why one 
must also not apply the pixel window correction to that 


same term. In fact, this statement is correct only for the 
shot noise, if the data are pixelized by the nearest-grid 
assignment (which we have adopted for our pixelizing 
scheme). This has been shown by Ref. [49] (see Eq. 20 of 
that work) for a 3-dimensional density field, but the same 
is true for a 2-dimensional field, as we are dealing with 
here. We have verified this using numerical simulations. 

In this paper, although we use maps at AT slde = 512 
(for which the maximum multipole is £ max = 1024), we 
restrict the analysis to Ct up to £ mazi ~ 500 where we 
have a reasonable signai-to-noise ratio. For these multi- 
poles the effect of the pixel window function is negligible, 
and thus we shall simplify our analysis pipeline by not ap- 
plying the pixel window correction to the observed power 
spectrum [60]. Therefore, our signal power spectrum es- 
timator is given by 

^signal _ CT7 /sky ~ 
t ^^bram^2 ’ ' ' 

where C N = (jV 7ipix )(l/A2 ix )/fi plx is the photon noise 
term, with A^pix, A p j x , and ft p i x the number of observed 
events, the exposure, and the solid angle, respectively, of 
each pixel, and the averaging is done over the unmasked 
pixels. We approximate the photon noise term by C n = 
(/) 2 47r/ s ky /-/V 7 , with N y denoting the total number of 
observed events outside the mask. This approximation 
is accurate at the percent level. Note that while (7| aw is 
always non- negative, it is possible for our estimator for 
the signal power spectrum C| lsnal to be negative due to 
the subtraction of the noise term. 

The beam window function in multipole space associ- 
ated with the full non-Gaussian PSF is given by 

WD m {E) = 2jt J 1 d cos 0P( (cos(0) )PSF (0 ; E), (5) 

where P£(cos(0)) are the Legendre polynomials and 
PSF(0; E) is the energy-dependent PSF for a given set 
of IRFs, with 0 denoting the angular distance in the dis- 
tribution function. The PSF used corresponds to the 
average for the actual pointing and live time history of 
the LAT and over the off-axis angle, as given by the gtpsf 
tool. We calculate the beam window functions for both 
the front- and back-converting events. 

The PSF of the LAT, and consequently the beam win- 
dow function, varies substantially over the energy range 
used in this analysis, and also non-negligibly within each 
energy bin. We treated this energy dependence by calcu- 
lating an average window function (W beam (£!;)) for each 
energy bin Ei , weighted by the intensity spectrum of the 
events in each bin, 

= -j- d E WD m {E) (6) 

'bin Ob 

where hit, s= d E {dN/dE) and E miUii and E maXii 

are the lower and upper edges of each energy bin. The 
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differential intensity dN/dE outside the mask in each 
map for the finely-gridded energy bins described in §11 
was used to approximate the energy spectrum for this 
calculation. 


D. Measurement uncertainties 

The la statistical uncertainty ac t on the measured an- 
gular power spectrum coefficients C'| lgnal is given by [50] 

__ / ^ / ^signal , @N 

Ct y (2^+l)/sky V 1 (VF beam ) 2 

where At is the width of the multipole bin (for binned 
data). 

After implementing the corrections for masking and for 
the beam window function to estimate the signal angular 
power spectrum via Eq. 4, the coefficients C| lgnal were 
binned in multipole with A^ = 50 and averaged in each 
multipole bin, weighted by the measurement uncertain- 
ties, 

(Ct) = (8) 

2^t l / a c l 

with Ci = C| lgn " 1 as calculated by Eq. 4 and ac t given 
by Eq. 7 with At = 1 and for 

the corresponding energy bin As expected, we find 
that the statistical measurement uncertainties calculated 
at the linear center of each multipole bin via Eq. 7 with 
At — 50 agree well with the scatter within each multipole 
bin. The value of Ct at multipoles 2 < t < 4 was found 
in most cases to be anomalously large [61], indicating 
the presence of strong correlations on very large angular 
scales, such as those that could be induced by the shape 
of the mask and by contamination from Galactic diffuse 
emission. To avoid biasing the value of the average C^ lgnal 
in the first bin by the values at these low multipoles, the 
multipole bins begin at t = 5. 

Finally, the angular power spectra of the front- and 
back-converting events were combined by weighted aver- 
aging, weighting by the measurement uncertainty on each 
data point. Due to the larger PSF associated with back- 
converting events, the measurement errors on the angular 
power spectra of the back-converting data set tend to be 
larger than those of the front-converting data set, par- 
ticularly at low energies and high multipoles where the 
suppression of the raw angular power due to the beam 
window function is much stronger for the back-converting 
data set. The difference between the measurement un- 
certainties associated with the front and back data sets 
is less prominent at higher energies. 



IV. EVENT-SHUFFLING TECHNIQUE 

One way to search for anisotropies is to first calculate 
the flux of particles from each direction in the sky (equal 


to the number of detected events from some direction 
divided by the exposure in the same direction), and then 
examine its directional distribution. The flux calculation, 
which requires knowledge of the exposure, depends on 
the effective area of the detector and the accumulated 
observation live time. 

The effective area, calculated from a Monte Carlo sim- 
ulation of the instrument, could suffer from systematic er- 
rors, such as miscalculations of the dependence of the ef- 
fective area on the instrument coordinates (off-axis angle 
and azimuthal angle). Naturally, any systematic errors 
involved in the calculation of the exposure will propagate 
to the flux, possibly affecting its directional distribution. 
If the magnitude of these systematic errors is compara- 
ble to or larger than the statistical power of the available 
data set, their effects on the angular distribution of the 
flux might masquerade as a real detectable anisotropy. 
For this reason, we cross-check our results using an alter- 
native method to construct an exposure map that does 
not rely on the Monte-Carlo-based calculation of the ex- 
posure implemented in the Science Tools. 
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FIG. 2: No-anisotropy sky map created by summing 20k 
shuffled maps using front- and back-converting events with 
E > 1 GeV, binned into a HEALPix map with N a ide = 256. 
The map projection is Hammer-Aitoff. The features in the 
no-anisotropy sky map result from the fact that the sky was 
not observed with uniform exposure. 

The starting point of this method is the construction 
of a sky map that shows how an isotropic sky would look 
as seen by the Fermi LAT. This sky map, hereafter called 
the “no-anisotropy sky map” , is directly proportional to 
the exposure map. 

One method of generating a no-anisotropy map is to 
randomize the reconstructed directions of the detected 
events (as in [41]). In the case that the angular distribu- 
tion of the flux is perfectly isotropic, a time-independent 
intensity should be detected when looking in any given 
detector direction. Possible time variation of the inten- 
sity would be due only to changes in the operating condi- 
tions of the instrument. A set of isotropic events can be 
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built by randomly coupling the times and the directions 
of real events in local instrument coordinates. The ran- 
domisation in this analysis was performed by exchanging 
the direction of a given real event in the LAT frame with 
the direction of another event selected randomly from the 
data set with uniform probability. Using this informa- 
tion, the sky direction is re-evaluated for the two events. 
By construction, the randomized data set preserves the 
exposure, the energy and angular (with respect to the 
LAT reference frame) distributions, and also accounts 
for the detector dead times. 

As already discussed in §11, for this analysis a cut of 
52° on the rocking angle was applied to limit possible 
photon contamination from the Earth’s albedo. For the 
shuffling technique, the analysis was performed with a re- 
duced field of view of the instrument, namely, the events 
used were selected to have an off-axis angle less than 50°. 
In this way, events with zenith angle exceeding 102° were 
removed. This selection cut avoids introducing asymme- 
tries in the exposure across the field of view due to cutting 
events based on zenith angle. 

The randomization was performed using the masked 
sky map described in §111, so that only real events with 
sky coordinates outside the masks were used, and the re- 
evaluated sky direction for each event was required to be 
in the unmasked region of the sky. This randomization 
process was repeated 20k times, separately for the front- 
and back-converting events, each time producing a shuf- 
fled sky map that is compatible with an isotropic source 
distribution. The final no-anisotropy sky map for each 
energy bin was produced by taking the average of these 
20k shuffled sky maps. For the available event statistics, 
averaging 20k shuffled maps was reasonably effective at 
reducing the Poisson noise associated with the average 
number of events per pixel. To reduce the number of re- 
quired shuffled maps by increasing the average number 
of events per pixel, the shuffled maps were constructed 
at slightly lower resolution (A s ide = 256) than was used 
in the default analysis. When analyzing the anisotropy 
with these exposure maps from the shuffling technique, 
count maps at Aside = 256 were used to construct the 
intensity maps. A no-anisotropy sky map is shown in 
Fig. 2. This sky map does not appear entirely uniform 
because the sky was not observed with uniform exposure. 

Although the no-anisotropy sky map is directly propor- 
tional to the exposure map, this method does not allow 
us to determine the absolute level of the exposure. We 
therefore constructed intensity maps (with arbitrary nor- 
malization) by dividing the real data maps in each energy 
bin by the no-anisotropj’ map for that energy bin, after 
first smoothing the no-anisotropy map with a Gaussian 
beam with a — 1° to reduce the pixel-to- pixel fluctuar 
tions due to the finite number of events available to use 
in the randomization. This smoothing beam size removes 
noise in the no-anisotropy sky map above £ ~ 200, and 
was chosen because we focus our search for anisotropies 
in that multipole range. Angular power spectra were 
then calculated from these intensity maps as in §111. Due 


to the arbitrary normalization of these intensity maps, 
we calculate only fluctuation angular power spectra of 
the data when using the exposure map produced by this 
shuffling technique. 

V. SIMULATED MODELS 

Detailed Monte Carlo simulations of Fermi LAT all- 
sky observations were performed to provide a reference 
against which to compare the results obtained for the real 
data set. The simulations were produced using the gto- 
bssim tool, which simulates observations with the LAT 
of an input source model. The gtobssim tool generates 
simulated photon events for an assumed spacecraft point- 
ing and live-time history, and a given set of IRFs. The 
P6.V3JDIFFUSE IRFs and the actual spacecraft point- 
ing and live-time history matching the observational time 
interval of the data were used to generate the simulated 
data sets. 

Two models of the gamma-ray sky were simulated. 
Each model is the sum of three components: 

1. GAL - a model of the Galactic diffuse emission 

2. CAT - the sources in the 1 1-month cata- 
log (1FGL) [46] 

3. ISO - an isotropic background 

Both models include the same CAT and ISO compo- 
nents, and differ only in the choice of the model for the 
GAL component. GAL describes both the spatial dis- 
tribution and the energy spectrum of the Galactic dif- 
fuse emission. The GAL component for the reference sky 
model used in this analysis (hereafter, MODEL) is the 
recommended Galactic diffuse model for LAT data anal- 
ysis, gll_iem-v02.f it [62], which has an angular resolu- 
tion of 0.5°. This model was used to obtain the 1FGL 
catalog; a detailed description can be found in Ref. [63]. 

An alternate sky model (ALT MODEL) was simulated 
for comparison, in order to test the possible impact of 
variations in the Galactic diffuse model. This model is 
internal to the LAT collaboration, and was built using 
the same method as gIlJ.em_v02 .f it, but differs pri- 
marily in the following ways: (i) this model was con- 
structed using 21 months of Fermi LAT observations, 
while gll_iem_v02.f it was based on 9 months of data; 
and (ii) additional large-scale structures, such as the 
Fermi bubbles [51], are included in the model through 
the use of simple templates. 

The sources in CAT were simulated with energy spec- 
tra approximated by single power laws, and with the 
locations, average integral fluxes, and photon spectral 
indices as reported in the 1FGL catalog. All 1451 
sources were included in the simulation. ISO repre- 
sents the sum of the Fermi-measured IGRB and an addi- 
tional isotropic component presumably due to unrejected 
charged particles; for this component the spectrum tem- 
plate isotropic-i enuv02.txt was used. 



10 


For both the MODEL and the ALT MODEL, the sum 
of the three simulated components results in a descrip- 
tion of the gamma-ray sky that closely approximates 
the angular-dependent intensity and energy spectrum of 
the all-sky emission measured by the Fermi LAT. Al- 
though the simulated models may not accurately repro- 
duce some large-scale structures, e.g., Loop I [52] and 
the Fermi bubbles, these features are not expected to in- 
duce anisotropies on the small angular scales on which 
we focus in this work. 


VI. RESULTS 

In this section we present the measured angular power 
spectra of the data, followed by the results of validation 
studies which examine the effect of variations in the de- 
fault analysis parameters, and by a comparison of the 
results for the data with those for simulated models. We 
summarize the main results of the angular power spec- 
trum measurements of the data and of key validation 
studies in Table II. 

Unless otherwise noted, the results are shown for data 
and models with the angular power spectra calculated 
after applying the default source mask which excludes 
sources in the 1FGL catalog and Galactic latitudes \b\ < 
30°. Due to the arbitrary normalization of the intensity 
maps calculated using the exposure map from shuffling, 
we show fluctuation angular power spectra for this data 
set. Intensity angular power spectra are presented for all 
other data sets. 

In the figures we show our signal angular power spec- 
trum estimator C| lgnal (Eq. 4), which represents the sig- 
nal after correcting for the power suppression due to 
masking, subtracting the photon noise, and correcting for 
the beam window function. A measurement that is in- 
consistent with zero thus indicates the presence of signal 
angular power. The (7| lgnaI shown is the weighted aver- 
age of this quantity for the maps of front and back events. 
The fluctuation angular power spectra C| lgnal /{/) 2 were 
calculated by dividing C| lgnal of the front and back events 
by their respective (/) 2 , and then averaging the angu- 
lar power spectra. For conciseness, in the figure labels 
Ct — G| aw // S k y is the raw angular power spectrum out- 
put by HEALPix corrected for the effects of masking. 
The error bars on points indicate the lcr statistical un- 
certainty in the measurement in each multipole bin as 
calculated by Eq. 7 with = 50 and with the bins be- 
ginning at l — 5. The binned data points are located at 
the linear center of each multipole bin. 


A. Angular power spectrum of the data 

We now present the results of the angular power spec- 
trum analysis of the data. We measure the angular power 
spectrum of the data after applying the default latitude 


cut and source mask, and refer to this as our default data 
analysis (DATA). We also measure the angular power 
spectrum of the data using the same masking and anal- 
ysis pipeline after performing Galactic foreground clean- 
ing, described below, and refer to this as the cleaned 
data analysis (DATA:CLEANED). These two measure- 
ments constitute our main results for the data, and so we 
discuss the energy dependence of the measured angular 
power (§VII) and present constraints on specific source 
populations (§VIII) for the results of both the default 
and cleaned data analyses. 

To minimize the impact of Galactic foregrounds we 
have employed a large latitude cut. However, Galactic 
diffuse emission extends to very high latitudes and may 
not exhibit a strong gradient with latitude, and it is thus 
important to investigate to what extent our data set may 
be contaminated by a residual Galactic contribution. For 
this purpose we attempt to reduce the Galactic diffuse 
contribution to the high-latitude emission by subtracting 
a model of the Galactic foregrounds from the data, and 
then calculating the angular power spectra of the residual 
maps. For the angular power spectrum analysis of the 
residual maps (cleaned data) we note that the noise term 
Cn is calculated from the original (uncleaned) map, since 
subtracting the model from the data does not reduce the 
photon noise level. 

In the following we use the recommended Galactic dif- 
fuse model gll_iem.v02.fit, which is also the default 
GAL model that we simulate, as described in §V. To tai- 
lor the model to the high-latitude sky regions considered 
in this work, the normalization of the model was adjusted 
by refitting the model to the data only in the regions 
outside the latitude mask. For the fit we used GaRDiAn 
which convolves the model with the instrument response 
(effective area and PSF). The normalization obtained 
in this way is, however, very close to the nominal one, 
within a few percent. 

We present the angular power spectra of the data be- 
fore and after Galactic foreground cleaning in Fig. 3; ex- 
panded versions of the angular power spectra for the 1- 
2 GeV and 2-5 GeV bins focusing on the high-multipole 
data are shown in Fig. 4. In both analyses, angular power 
at l > 155 is measured in the data in all energy bins con- 
sidered, and the angular power spectra for the default 
and cleaned data are in good agreement in this multipole 
range. In the default data, the large increase in angular 
power at l < 155 in the two energy bins spanning 1- 
5 GeV is likely due to contamination from the Galactic 
diffuse emission which features correlations on large an- 
gular scales, but may also be attributable in part to the 
effects of the source mask (see §V1F). 

At i > 155 the measured angular power does not ex- 
hibit a clear scale dependence in any energy bin. The 
results of fitting the unbinned signal angular power spec- 
trum estimator for 155 < i < 504 in each energy bin to 
a power law C| lgnal a (£/<o) n with £q = 155 are given 
in Table I for the default data analysis. In each energy 
bin, the angular power spectrum for 155 < i < 504 is 
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FIG. 3: Comparison of intensity angular power spectra of the data and Galactic-foreground-cleaned data. For £ > 155 the 
measured power at all energies is approximately constant in muitipole, suggesting that it originates from one or more unclustered 
source populations. The large increase in angular power in the default data at £ < 155 in the 1-2 and 2-5 GeV bins is likely 
attributable largely to contamination from Galactic diffuse emission. In these two energy bins, foreground cleaning primarily 
reduces angular power at £ < 155, with the most significant reductions at £ < 105. At energies greater than 5 GeV the effect 
of foreground cleaning is small for £ > 55. Expanded versions of the top panels are shown in Fig. 4. 


consistent with a Poisson spectrum (constant in multi- 
pole, i.e., n = 0 falls within the 95% CL range of the 
best-fit power-law index), as expected for the angular 
power spectrum of one or more unclustered source pop- 
ulations. However, we emphasize that the uncertainty in 
the scale dependence is appreciable, particularly for the 
10-50 GeV bin. 

TABLE I: Multipole dependence of intensity angular power in 
the data (default analysis) for 155 < £ < 504 in each energy 
bin. The best-fit power-law index n in each energy bin is 
given with the associated \ 2 per degree of freedom (d.o.f.) of 
the fit. 


-Emin 

Em&x 

n 

X 2 /d.o.f. 

1.04 

1.99 

-1.33 ±0.78 

0.38 

1.99 

5.00 

-0.07 ± 0.45 

0.43 

5.00 

10.4 

-0.79 ± 0.76 

0.37 

10.4 

50.0 

-1.54 ±1.15 

0.39 


In light of the scale independence of the angular power 


at £ > 155, we associate the signal in this multipole range 
with a Poisson angular power spectrum and determine 
the best-fit constant value of the angular power Cp and 
the fluctuation angular power Cp/(/) 2 over 155 < £ < 
504 in each energy bin, by weighted averaging of the un- 
binned measurements. These results for the default and 
cleaned data are summarized in Table II, along with the 
results obtained for the data using an updated source cat- 
alog to define the source mask and for a simulated model, 
which will be discussed in §VIG and §VIH, respectively. 

We note that the associated measurement uncertain- 
ties can be taken to be Gaussian, in which case the re- 
ported significance quantifies the probability of the mea- 
sured angular power to have resulted by chance from a 
truly uniform background. We consider a 3a or greater 
detection of angular power (Cp) in a single energy bin 
to be statistically significant. For the default data, the 
best-fit values of Cp indicate significant detections of an- 
gular power in the 1-2, 2-5, and 5-10 GeV bins (6.5(7, 
7.2(7, and 4.1(7, respectively), while in the 10-50 GeV bin 
the best-fit Cp represents a 2.7cr measurement of angu- 
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FIG. 4: Expanded versions of top panels of Fig. 3, focusing on the high-multipole angular power. 


TABLE II: Best- fit values of the angular power Cp and fluctuation angular power Cp/{/) 2 in each energy bin over the multipole 
range 155 < £ < 504. Results are shown for the data processed with the default analysis pipeline, the foreground-cleaned data, 
the data analyzed with the 2FGL source mask, and the default simulated model. Significance indicates the measured angular 
power expressed in units of the measurement uncertainty a; the measurement uncertainties can be taken to be Gaussian. 



jFmin 

[GeV] 

flmax 

[GeV] 

C P 

[(cam -2 s _ 1 sr" 1 ) 2 sr] 

Significance 

Cp/(I ) 2 

[10 6 sr] 

Significance 

DATA 

1.04 

1.99 

7.39 ±1.14 x 10” 18 

6.5 a 

10.2 ± 1.6 

6.5a 


1.99 

5.00 

1.57 ±0.22 x 10" 18 

7.2 u 

8.35 ± 1.17 

7.1(7 


5.00 

10.4 

1.06 ± 0.26 x 10" 19 

4.1(7 

9.83 ± 2.42 

4.1(7 


10.4 

50.0 

2.44 ±0.92 x 10" 20 

2.7(7 

8.00 ± 3.37 

2.4(7 

DATA: CLEANED 

1.04 

1.99 

4.62 ±1.11 x 10" 18 

4.2<7 

6.38 ± 1.53 

4.2(7 


1.99 

5.00 

1.30 ±0.22 x 10" 18 

6.0<7 

6.90 ± 1.16 

5.9a 


5.00 

10.4 

8.45 ±2.46 x 10" 20 

3.4(7 

8.37 ± 2.41 

3.5(7 


10.4 

50.0 

2.11 ±0.86 x 10" 2 ° 

2.4<7 

7.27 ± 3.36 

2.2(7 

DATA:2FGL 

1.04 

1.99 

5.18 ± 1.17 x 10" 18 

4.4<7 

7.23 ± 1.61 

4.5(7 


1.99 

5.00 

1.21 ±0.28 x 10" 18 

5.3a 

6.49 ± 1.22 

5.3(7 


5.00 

10.4 

8.38 ± 2.72 x 10“' 20 

3.1(7 

7.67 ± 2.54 

3.0cr 


10.4 

50.0 

8.00 ± 9.57 x 10" 21 

0.8(7 

2.28 ± 3.52 

0.6<7 

MODEL 

1.04 

1.99 

1.89 ±1.08 x 10" 18 

0.7a 

2.53 ± 1.47 

1.7(7 


1.99 

5.00 

1.92 ± 2.10 x 10" 19 

0.9(7 

0.99 ± 1.12 

0.9a 


5.00 

10.4 

3.41 ± 2.60 x 10" 20 

1.3(7 

3.04 ± 2.34 

1.3(7 


10.4 

50.0 

0.62 ± 9.63 x 10" 21 

0.1(7 

0.24 ± 3.02 

0.1(7 


lar power. We further note that the best-fit value of the 
fluctuation angular power over all four energy bins (see 
§VII and Table IV) yields a detection with greater than 
lOcr significance for the default data. 

For the 1-2 GeV and 2-5 GeV energy bands the clean- 
ing procedure results in a significant decrease in the an- 
gular power at low multipoles (£ < 105), and a smaller 
reduction at higher multipoles. However, the decrease is 
small for £ > 155, and angular power is still measured at 
all energies, at slightly lower significances (see Table II). 
We emphasize that the detections in the three energy 
bins spanning 1-10 GeV remain statistically significant, 
and the best-fit fluctuation angular power over all en- 


ergy bins is detected at greater than 8<r significance. For 
energies above 5 GeV the foreground cleaning does not 
strongly affect the measured angular power spectrum for 
t > 55. At all energies the decrease in angular power 
at low multipoles can be attributed to the reduction of 
Galactic foregrounds which feature strong correlations on 
large angular scales. We conclude that contamination of 
the data by Galactic diffuse emission does not have a sub- 
stantial impact on our results at the multipoles of inter- 
est (£ > 155). This conclusion is in agreement with that 
of Ref. [39], which found that the Galactic foregrounds 
have a rapidly declining angular power spectrum above 
£- 100 . 
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To further study the expected angular power spectrum 
of Galactic foregrounds, we analyzed the angular power 
spectrum of the E(B-V) emission map of Ref. [53] (here- 
after SFD map), which is proportional to the column 
density of the interstellar dust, after masking |b| < 30 c 
as in our default analysis. The SFD map is a good tracer 
of the Galactic interstellar medium (ISM) away from the 
Galactic plane, the spatial structure of which should be 
reflected in the diffuse gamma-ray emission produced by 
interactions of cosmic rays with the ISM. It has an an- 
gular resolution of 6 arcminutes, much smaller than the 
intrinsic resolution of the GAL model map (^0.5°), and 
smaller than the map resolution used in this study, and 
so it accurately represents the small-scale structure of 
the ISM on the angular scales accessible to this analy- 
sis. We found that the SFD map produces an angular 
power spectrum with a slightly harder slope than the 
default GAL model, and consequently features more an- 
gular power at high multipoles. However, like the GAL 
model the SFD map angular power spectrum falls off 
quickly with multipole compared to a Poisson spectrum, 
and the amplitude of the SFD map angular power is be- 
low that measured in the data for t > 100. This further 
reinforces the conclusion that Galactic foreground con- 
tamination cannot explain the observed high-multipole 
angular power in the data. 


B. Validation with a simulated point source 
population 

To ensure that our analysis procedure accurately re- 
covers an input angular power spectrum, and in particu- 
lar that the result is not biased by instrumental effects, 
we compare the angular power spectrum calculated for 
a simulated point source population with the theoreti- 
cal prediction for that population. It is straightforward 
to calculate the expected angular power spectrum of un- 
clustered point sources, once a flux distribution function, 
dN/dS (in units of cm 2 s sr _1 ), and a source detec- 
tion flux threshold, S c (in units of cm -2 s"* 1 ), are pro- 
vided. The angular power spectrum of an unclustered 
point source population is the Poisson component of the 
angular power Cp, which takes the same value at all mul- 
tipoles and is given by 

<9) 

For our source population model we adopt the best- 
fit flux distribution for the high-latitude Fermi sources, 
reported in [10], which describes dN/dS with a broken 
power-law model: 

^ = AS^, S > Sb, 

= AS; 0l+l) *S-h, S < S b , (10) 

which contains four free parameters, A, 0i, 0 2} and S b . 
For this form of dN/dS, the source power spectrum can 



FIG. 5: Predicted amplitude of the source angular power 
spectrum, Cp (see Eq. 9), for energies of 1.04-10.4 GeV as 
a function of a source detection threshold flux, S c • 



FIG. 6: Intensity angular power spectrum of a simulated ob- 
servation of the source population model, compared with the 
theoretical prediction ( shaded band). The angular power spec- 
trum of the simulated population is in excellent agreement 
with the prediction over a large multipole range. 


be found analytically (for S c > S 5): 



— 1 

Xo 

1 

Xo 

to 

1 

?6\ 3 ' /3l l 

3- A V 

U 


(id 


A fit for the simulated source population for 1.04- 
10.4 GeV yields A = (1.90 ± 0.48) x 10" 13 (180/t r) 2 , 
0 1 = 2.213 ± 0.073, 02 = 1 533 ± 0.007, and S b = 
1.41 x 10” 9 cm -2 s -1 . Note that the errors are cor- 
related. Fig. 5 shows the predicted Cp for this source 
model for threshold fluxes in the range S c = 0.8 x 10“ 7 - 
4.2 x 10 -7 cm -2 s -1 . The error bars are calculated from 
the full covariance matrix of the above parameters. Al- 
though we have used zero as the lower limit of the integral 
in Eq. 9, using the actual lower limit of the flux distri- 
bution adopted for the simulated population results in a 
negligible difference in the predicted Cp. 
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We simulated this source population model with gtob - 
ssim using the same procedure as described in §V. The 
simulated population comprises nearly 20k point sources 
distributed randomly across the entire sky, with each 
source’s flux drawn from the flux distribution specified 
above. The photon spectrum of each source is modeled 
as a power law with a spectral index T (dN/dE oc E~ Tg ) 
drawn from a Gaussian distribution with mean of 2.40 
and a standard deviation of 0.28. The simulated events 
were processed and the angular power spectrum of this 
source model calculated using the same procedure as was 
used for the data and other simulations in this study, ex- 
cept that the energy range of the map was chosen to be 
1.04-10.4 GeV, and no mask was applied. 

The fluxes of the ~ 20k simulated sources were drawn 
from a flux distribution in which the maximum possi- 
ble flux (E > 100 MeV) that could be assigned to a 
source was 10 _5 cm -2 s -1 , however the maximum flux 
of any source in the simulation, which represents a sin- 
gle realization of this source population, was ^3x 10^ 6 
cm -2 s -1 . We take these values as the upper and lower 
bound on the source detection threshold flux (E > 100 
MeV) corresponding to the simulated model, since we 
do not impose a source detection threshold by masking 
or otherwise excluding simulated sources above a specific 
threshold flux. A spectral index V = 2.4 is assumed to 
determine the threshold fluxes in the 1.04-10.4 GeV en- 
ergy band. Prom these threshold fluxes we calculate the 
corresponding upper and lower bound on the predicted 
Cp in the 1.04-10.4 GeV energy band. 

The angular power spectrum for the simulated source 
population calculated via the analysis pipeline used in 
this study is presented in Fig. 6, with the shaded region 
indicating the predicted range of Cp (the mean values of 
Cp at the upper and lower flux threshold); for a given 
model Cp is independent of multipole, thus we expect 
the recovered angular power spectrum to be independent 
of multipole with amplitude within the shaded region. 
The angular power spectrum recovered from the simu- 
lated data is in excellent agreement with the prediction 
up to multipoles of £ ~ 800. Above £ ~ 800, the upturn 
in the measured angular power spectrum is likely due to 
inaccuracies in the modeling of the beam window func- 
tion, which can introduce features on very small angular 
scales. In the remainder of this study, we present results 
only for the multipole range £ = 5 to £ = 504. 


C. Sensitivity to the exposure map calculation 

To investigate the possibility that potential inaccura- 
cies in the exposure map calculation for the default anal- 
ysis might generate spurious anisotropy in the intensity 
maps, we compare the fluctuation angular power spec- 
tra of the data using our default analysis pipeline with 
the results obtained after replacing the default exposure 
map with that generated by the event shuffling technique 
described in §IV. This is shown in Figs. 7 and 8. In 


these two figures only, the results from the default data 
analysis were obtained from maps of HEALPix resolution 
iV g id e = 256 to match the resolution of the maps using 
the exposure determined from the shuffling technique. 
All other results presented in this study were obtained 
from N s ide = 512 maps. Due to the reduced map resolu- 
tion, the pixel window function has a small effect on the 
angular power spectra shown in Figs. 7 and 8, however 
it affects the results of the default analysis and the anal- 
ysis using the shuffled exposure map in the same way, 
and so these results can still be directly compared for 
the purpose of checking the effect of the exposure map 
calculation. 

The results of the two analysis methods are in good 
agreement at all energies and multipoles considered, ex- 
cept for slight deviations at £ < 55 for 1-5 GeV. We 
caution that at these low multipoles the measured angu- 
lar power spectra may be strongly affected by the mask, 
which has features on large angular scales. The slight 
differences in the data selection cuts for the analysis us- 
ing the exposure map from the shuffling technique com- 
pared to those for the default data analysis could lead 
to the observed differences in the low-multipole angular 
power spectra. The differences could also result from sys- 
tematics in the Monte-Carlo-based exposure calculation 
implemented in the Science. Tools, leading to inaccura- 
cies in the exposure map which vary on large angular 
scales. As we do not focus on the low-multipole an- 
gular power in this study, we defer a full investigation 
of this issue to future work. The agreement at £ > 55 
demonstrates that any potential spatially-dependent in- 
accuracies in the Science Tools exposure calculation have 
a negligible impact on the angular power spectra in the 
multipole range of interest. In particular, from the con- 
sistency of the two methods we conclude that using the 
Monte-Carlo-based exposure calculation does not induce 
spurious signal anisotropy in our results. 


D. Dependence on the PSF model 

We examine the impact of variations in the assumed 
PSF on the results of the analysis by comparing the beam 
window functions (Eq. 5) for the PSF implemented in the 
P6JV3 IRFs used in this analysis to those for the PSF in 
the more recently updated P6-V11 IRFs. The P6-V11 
IRFs use a modified functional form for the PSF, and for 
energies above 1 GeV the PSF implemented in P6-V11 
was calibrated using in-flight data, while in P6-V3 the 
PSF was based on Monte Carlo simulations. Fig. 9 shows 
the beam window functions for the PSF associated with 
the front- and back-converting events for each set of IRFs, 
at the log center of each energy bin used in this analy- 
sis. The small variation between the window functions 
of the two IRFs confirms that differences between the 
PSF models in these two IRFs are not large enough to 
affect the anisotropy measurement on the angular scales 
to which this analysis is sensitive. 
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FIG. 7: Fluctuation angular power spectra Ce/{I) 2 calculated using the default analysis pipeline compared with those obtained 
using the exposure map from the event shuffling technique described in §TV. Angular power is measured in all four energy bins 
by botn analysis methods. The lack of significant differences at the multipoles of interest between the angular power spectra 
yielded by the two methods demonstrates that any inaccuracies in the exposure map have a negligible impact on the measured 
angular power spectra. Expanded versions of the top panels are shown in Fig. 8. 
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FIG. 8: Expanded versions of top panels of Fig. 7. 


E. Dependence on the latitude mask 


In this analysis we apply a generous latitude mask 
to reduce contamination of the data by Galactic diffuse 


emission. The mask is intended to remove enough con- 
tamination so that the measured angular power can be 
attributed to sources whose distribution is statistically 
isotropic in the sky region we consider, i.e., a distribution 
which does not show any preferred direction on the sky. 
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Multipole / 



FIG. 9: Comparison of the beam window functions for the P6_V3 and P6-V11 IRFs; the P6-V3 IRFs are the default used in this 
analysis. The quantity Wf , which is the factor by which the angular power is suppressed due to the finite angular resolution 
of the instrument, is shown for the front-converting ( left panel) and back-converting ( right panel) events, evaluated at the 
log-center of each energy bin used in this analysis. The differences between the Wl of these two IRFs are small (< few percent) 
at all energies considered, indicating that our results are insensitive to the differences between the PSF models implemented in 
these IRFs. 


In particular, we wish to exclude sources whose angular 
distribution exhibits a strong gradient with Galactic lat- 
itude. The effectiveness of the mask at reducing the con- 
tribution to the angular power from a strongly latitude- 
dependent component can be evaluated by considering 
the angular power spectrum of the data as a function of 
latitude cut. The results are shown in Figs. 10 and 11. 

At low multipoles (£ < 100), increasing the latitude cut 
significantly reduces the angular power, indicating that 
in this multipole range the contamination by a strongly 
latitude-dependent component, such as Galactic diffuse 
emission, is considerable. For 155 < t < 254 at 1-2 GeV 
and 2-5 GeV, the angular power measured using the 30° 
latitude mask is noticeably smaller than when using the 
20° latitude mask. However, at all energies there are 
no significant differences in the angular power measured 
for t > 155 using the 30° and 40° latitude masks, and 
for energies greater than 5 GeV the 20° latitude mask 
also yields consistent results. We conclude that apply- 
ing the 30° latitude mask is sufficient to ensure that 
no significant amount of the measured angular power at 
£ > 155 originates from the Galactic diffuse emission or 
from any source class that varies greatly in the region 
30° < |6| < 40°. 


F. Effects of masking on the power spectrum 

To verify that the results do not depend sensitively 
on the angular radius of the source mask, in Figs. 12 
and 13 we compare the results when masking a 1° angular 
radius around each source with those when masking the 
2° radius used as the default in this work. 

In the 1-2 GeV energy bin the results show significant 
differences at £ < 155, however for £ > 155 (the multipole 


range of interest) the angular power spectra for the 1° 
and 2° source mask cases agree within the error bars. In 
the higher energy bins the angular power spectra in all 
except the first multipole bin (5 < £ < 55, well below 
the range of interest) agree within the error bars. Since 
varying the angular size of the region masked around 
each source does not significantly change the measured 
angular power at £ > 155, we conclude that any features 
that may be induced in the angular power spectra by 
the morphology of the source mask are confined to low 
multipoles and therefore do not affect the measurements 
of Cp reported in this work. 

In addition, we have confirmed that the angular power 
spectra of the front- and back-converting events are in 
good agreement within each energy bin in the multipole 
range of interest (£ > 155), and are generally consistent 
at £ < 155 even in the 1-2 GeV energy bin where the 95% 
containment radius of the PSF of the back-converting 
events is comparable to the angular radius used for the 
source mask. Consequently, although the PSF associated 
with the back-converting events is larger than that of the 
front-converting events, the consistency of their angular 
power spectra implies that the source masking is suffi- 
ciently effective even at low energies. 

The sharp latitude cut used in this analysis also has the 
potential to induce features in the angular power spec- 
trum, although these would be expected to appear on the 
large angular scales characteristic of the morphology of 
the mask. We therefore note that the stability of the an- 
gular power spectra at £ > 155 for latitude cuts masking 
at least |6| < 30°, discussed in §VIE and demonstrated 
in Figs. 10 and 11, indicates that the latitude mask does 
not induce features in the power spectrum at the angular 
scales of interest. 

The analysis of the simulated isotropic component, 
presented in §VIH and in Figs. 18 and 19, provides an- 
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FIG. 10: Intensity angular power spectra of the data calculated with different latitude cuts. The point source mask was applied 
in addition to the latitude mask in all cases. The differences between the results masking |6| < 30° (the default latitude cut) 
and |6] < 40 c are small for £ > 155 for all four energy bins, demonstrating that the power observed in the data at these 
multipoles is not strongly correlated with a component that has a strong latitude dependence in the range 30° < |6| < 40°, 
such as the Galactic diffuse emission. At energies above 5 GeV convergence is seen for multipoles £ > 155 even when masking 
only |6| < 20°. Points from different data sets are offset slightly in multipole for clarity. The lowest multipole data point for 
the |ft| < 20° mask in each panel is above the range shown in the figure. Expanded versions of the top panels are shown in 
Fig. 11. 




FIG. 11: Expanded versions of top panels of Fig. 10. 
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FIG. 12: Intensity angular power spectra of the data calculated with a mask excluding a 1° or 2° angular radius around each 
source; excluding a 2° angular radius is the default in this analysis. The default latitude mask excluding |b| < 30° was applied 
in addition to the source mask in all cases. At all energies the angular power spectra obtained using the different source mask 
radii are consistent at £ > 155 (the multipole range of interest), and above 2 GeV the results are consistent at t > 55. Expanded 
versions of the top panels are shown in Fig. 13. 
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FIG. 13: Expanded versions of top panels of Fig. 12. 


other means of assessing the impact of the mask on the 
angular power spectra. Since the isotropic component 
should only contribute to the monopole (£ = 0) term of 
the power spectrum, statistically significant deviations 
from zero power at i > 0 can be attributed to the use 


of the mask. We emphasize that the consistency of the 
angular power of the isotropic component with zero at 
i > 155 indicates that, despite the complex morphology 
of the total mask, the mask does not induce features in 
the angular power spectrum at the multipoles of interest 
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[t > 155). 


G. Dependence on the set of masked sources 

The recently-released second Fermi LAT source cata- 
log (2FGL) [54] is an update to the 1FGL catalog used 
to define the default source mask adopted in this v/ork. 
The 2FGL catalog reports the detection of 1873 sources, 
compared to the 1451 included in the 1FGL catalog. 

We briefly comment that one motivation for using the 
1FGL catalog, rather than the 2FGL catalog, to define 
the source mask in our default analysis is that the 1FGL 
catalog was also used in the Fermi LAT source count 
distribution analysis [10]. The results of that study are 
closely related to the interpretation of the results of the 
current analysis, and so our choice to mask that same 
source list in our default analysis allows the results of 
the two analyses to be used together straightforwardly. 
However, it is natural to ask to what extent the measured 
angular power reported in the data may be attributable 
to the additional sources resolved in the 2FGL catalog. 

We address this question by analyzing the data using 
a source mask defined by the 2FGL sources and compar- 
ing the results to those obtained using the 1FGL source 
mask. We repeat the analysis of the data using the de- 
fault pipeline, changing only the source mask; the total 
mask is defined by the source mask combined with the 
default latitude cut masking |5| < 30°. When combined 
with the default latitude cut, the 2FGL source mask re- 
sults in an unmasked sky fraction f s y y = 0.295, a small 
decrease compared to f & ky = 0.325 when using the 1FGL 
source mask. 

The angular power spectra of the data analyzed using 
the 2FGL catalog to define the source mask are shown 
in Figs. 14 and 15, compared with the results of the de- 
fault data analysis which uses the 1FGL catalog. The 
angular power Cp measured in the data using the 2FGL 
source mask is reduced relative to the 1FGL case (see 
Table II), while the measurement uncertainties remain 
roughly the same as in the 1FGL case. The decrease in 
Cp is ~ 20-30% in the 1-2, 2-5, and 5-10 GeV energy 
bins, however significant detections (> 3 a) are still found 
in these three bins. A ~ 70% decrease in Cp is seen in 
the 10-50 GeV bin, and due to the large measurement 
uncertainty the significance of the measurement in this 
bin falls from 2.7a to 0.8a. The significance of the de- 
tected fluctuation angular power over all four energy bins 
remains greater than 7a. 

We can estimate the expected decrease in angular 
power when masking the 2FGL sources by calculating the 
difference in angular power produced when the source de- 
tection threshold is reduced from the 1FGL to the 2FGL 
catalog level, following the approach used to calculate 
the angular power of the simulated sources in §VIB. 
We assume the sources follow the flux distribution func- 
tion dN/dS given by Eq. 10 with the same parameters 
given in that section as the best-fit for the high-latitude 


Fermi sources in the 1.04-10.4 GeV band. Calculating 
Cp via Eq. 11 for an assumed flux threshold of ~ 5 x 
10“ lo cm“ 2 s” 1 , appropriate for the 1FGL catalog [46], 
yields Cp ~ 9.4 x 10“ 18 (cm' 2 s _1 sr -1 GeV -1 ) 2 sr. Us- 
ing a lower flux threshold of ^ 4 x 10 _10 cm -2 s _1 , appro- 
priate for the 2FGL catalog [54], gives Cp ~ 6.8 x 10 -18 
(cm -2 s _1 sr -1 GeV -1 ) 2 sr, which is indeed a roughly 
30% decrease in Cp, as observed in the data. 


H. Comparison of data and simulated models 

To understand the origin of the angular power mea- 
sured in the data, we compare the angular power spec- 
tra of the default data to those of the default simulated 
model and the alternate simulated model, described in 
§V. The simulated models were processed and their an- 
gular power spectra calculated using the same analysis 
pipeline as the data, and thus we expect the angular 
power spectra of the data and models to be consistent 
if the models accurately reflected the statistical proper- 
ties of the emission on the relevant angular scales. 

Figs. 16 and 17 present the angular power spectra of 
the data and models. The angular power spectra of the 
two models agree very well at all energies at multipoles 
above t = 105. At all energies and scales, both models 
exhibit less angular power than the data. Moreover, the 
amplitude of the detected angular power in both models 
is inconsistent with that of the data at > 95% CL in the 
three energy bins spanning 1-10 GeV, and at > 90% CL 
in the 10-50 GeV bin (see Table III). The lack of signifi- 
cant power at high multipoles in either simulated model 
indicates that the Galactic diffuse emission, as imple- 
mented in these models, provides a negligible contribu- 
tion to the anisotropy i > 155. At lower multipoles, the 
discrepancy between the data and models and between 
the two models may be due to the presence of large-scale 
features in the data which are not included in the mod- 
els, however we defer a full investigation of the origin of 
the low- multipole angular power to future v/ork. 

TABLE III: Significance of the difference ACp between in- 
tensity angular power Cp for 155 < t < 504 in the default 
data and the default simulated model in each energy bin. 
The associated measurement uncertainties can be taken to 
be Caussian. 


Emin 

E m ax 

Significance of ACp 

1.04 

1.99 

3.5<r 

1.99 

5.00 

4.5a 

5.00 

10.4 

2.0a 

10.4 

50.0 

1.7a 


The contributions to the angular power spectrum of 
the individual components of the default model are shown 
in Figs. 18 and 19. At all energies the only component 
contributing significantly to the total power is the Galac- 
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FIG. 14: Intensity angular power spectra of the data calculated using the source mask defined by the 2FGL catalog compared 
with the results using the 1FGL catalog; the source mask defined by the 1FGL catalog is the default used in this analysis. 
The angular power at £ > 155 is smaller in the 2FGL case by ~ 20-30% in the bins spanning 1-10 GeV and by ~ 70% at 
10-50 GeV. Expanded versions of the top panels are shown in Fig. 15. 
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FIG. 15: Expanded versions of top panels of Fig. 14. 


tic diffuse emission. The contribution from the isotropic 
component is negligible, since this component is isotropic 
by construction and thus, after the photon noise is sub- 
tracted, it should only contribute to the monopole {t = 0) 
term. The deviations from zero of the isotropic compo- 


nent in the lowest multipole bin (5 < t < 54) may be due 
to imperfect correction of the effects of the mask in this 
multipole regime. The source catalog component con- 
tributes zero power at all energies and multipoles since 
the emission maps of this simulated component contain 
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FIG. 15: Angular power spectra of the data, the default simulated model (MODEL), and the alternate simulated model (ALT 
MODEL). The angular power spectra of the two models are in good agreement in all energy bins. The smaller amplitude 
angular power at i > 155 measured at lower significance in both models is inconsistent with the angular power observed in the 
data a: all energies. Points from different data sets are offset slightly in multipole for clarity. Expanded versions of the top 
panels axe shown in Fig. 17. 
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FIG. 17: Expanded versions of top panels of Fig. 16. 


onlv events from sources which are masked in the analy- 
sis. The consistency of the source catalog angular power 
with zero indicates that the source masking is effective. 
We remark that in general the angular power spectra 
of distinct components are not linearly additive due to 


contributions from cross-correlations between the com- 
ponents. The total power of the model is, however, very 
consistent with the total power in the Galactic compo- 
nent, with slight discrepancies likely arising from mask- 
ing effects, since the Galactic and isotropic components 
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FIG. 13: Angular power spectra of the components of the default simulated model (MODEL). As expected, most of the total 
angular power at all multipoles (TOTAL MODEL) is due to the GAL component. By construction, the isotropic component 
(ISO) component contributes no significant angular power; likewise, the source component (CAT) provides no contribution 
because all sources were masked. Points corresponding to the TOTAL MODEL are offset slightly in multipole for clarity. 
Expanded versions of the top panels are shown in Fig. 19. 




FIG. 19: Expanded versions of top panels of Fig. 18. 


should have no cross-correlation power and the simulated VII. ENERGY DEPENDENCE OF 

sources were fully masked. ANISOTROPY IN THE DATA 


The energy dependence of the fluctuation angular 
power can be used to identify the presence of multiple 


23 


TABLE IV: Energy dependence of angular power for 155 < £ < 504 in each energy bin for the data processed with the 
default analysis pipeline and the Galactic-foreground-cleaned data. The best-fit constant value of the fluctuation angular 
power (Cp/{/) 2 ) over 1-50 GeV is obtained by weighted averaging of Cp/(I) 2 of the four energy bins. The best-fit parameters 
and associated x 2 P er degree of freedom (d.o.f.) are given for fits of the fluctuation angular power to Cp/{I ) 2 = Af(E/Eo)~ Fp 
and the differential intensity angular power to Cp/(A£) 2 = Ai(E/Eo)~ Vl : with Eq = 1 GeV. The value of A\ is given in terms 
of Ai/Ai.o where Ai,o = 10 -18 (cm -2 s" 1 sr -1 GeV -1 ) 2 sr. 



(Cp /( if) 

[10~ 6 sr] 

Af 

[10- 6 sr] 

r F 

X a /d.o.f. 

Ai/Ai.o 

Ti 

X 2 /d.o.f. 

DATA 

9.05 ± 0.84 

9.85 ± 1.73 

0.076 ± 0.139 

0.41 

45.1 ± 7.8 

4.79 ±0.13 

0.19 

DATA : CLEANED 

6.94 ± 0.84 

6.31 ± 1.44 

-0.082 ± 0.158 

0.12 

29.4 ± 6.6 

4.66 ± 0.15 

0.035 
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FIG. 20: Anisotropy energy spectra of the data. Top : Fluc- 
tuation anisotropy energy spectrum. The data are consistent 
with no energy dependence over the energy range considered, 
although a mild energy dependence is not excluded. Bottom : 
Differential intensity anisotropy energy spectrum. The energy 
dependence is consistent with that arising from a single source 
population with a power-law intensity energy spectrum with 
spectral index r 8 = 2.40±0.07 for the default data (2.33±0.08 
for the cleaned data). 


distinct contributors to the emission [45]. Because the 
fluctuation angular power characterizes only the angular 
distribution of the emission, independent of the intensity 
normalization, it is exactly energy- independent for a sin- 
gle source class as long as the members of the class have 


the same observed energy spectrum. In general, the fluc- 
tuation angular power of a single source class may show 
energy dependence due to large variation of the energy 
spectra of individual sources within a population, and, 
for cosmological source classes, the effects of redshifting 
and attenuation of high-energy gamma rays by the extra- 
galactic background light (EBL). Redshifting and EBL 
attenuation is expected to be important only for popula- 
tions for which a significant fraction of the observed in- 
tensity originates from high-redshift members, with EBL 
attenuation relevant only at observed energies of several 
tens of GeV. All of these effects are most prominent when 
the source spectra have hard features such as lines or cut- 
offs; smoothly- varying source spectra, such as power-law 
energy spectra, typically generate more mild energy de- 
pendence in the fluctuation angular power. 

The fluctuation anisotropy energy spectrum of the 
data is shown in the top panel of Fig. 20. The fluctuation 
angular power Cp/(/) 2 in each energy bin was obtained 
by weighted averaging of the unbinned fluctuation angu- 
lar power spectrum over 155 < i < 504, weighting the 
measured angular power at each multipole by its mea- 
surement uncertainty; these values are reported in Ta- 
ble II. Each point is located at the logarithmic center of 
the energy bin. 

A power-law fit of the fluctuation angular power as a 
function of energy Cp/(/) 2 oc E~ Tf yields Tp = 0.076 ± 
0.139 (-0.082 ± 0.158 for the cleaned data), consistent 
with no energy-dependence over the energy range con- 
sidered. The best-fit constant value of Cp/(/) 2 across all 
four energy bins is 9.05±0.84 x 10 -6 sr (6.94±0.84x 10 -6 
sr for the cleaned data). The results of these fits for the 
data with and without foreground cleaning are summa- 
rized in Table IV, along with the results for the energy de- 
pendence of the intensity angular power, discussed below. 
The. lack of a clear energy dependence in the fluctuation 
angular power is consistent with a single source class pro- 
viding the dominant contribution to the anisotropy and 
a constant fractional contribution to the intensity over 
the energy range considered, although due to the large 
measurement uncertainties contributions from additional 
source classes cannot be excluded. This is especially true 
for sources whose contribution to the intensity peaks at 
E > 10 GeV. Furthermore, due to the coarseness of the 
energy binning, this analysis is not sensitive to features 
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in the anisotropy energy spectrum localized to narrow 
energy bands. 


If a single source class dominates the anisotropy at 
all energies considered, the differential intensity angu- 
lar power spectrum <7// (A#) 2 scales with energy as the 
intensity energy spectrum squared (dN/dE ) 2 of that 
source class. For example, for a source class with a power- 
law photon spectrum dN/dE a E~ r \ Ct/(AE ) 2 a 
E~ 2r *. We can therefore use this energy scaling to con- 
strain the energy spectrum of the dominant contributor 
to the anisotrop 3 r , under the assumption that the mea- 
sured angular power (but not necessarily the total mea- 
sured intensity) originates from a single source class. 


Here we obtain the differential intensity angular power 
Cp/(AE ) 2 by dividing the intensity angular power Cp in 
each energy bin by the bin size squared. The differential 
intensity anisotropy energy spectrum of the data is shown 
in the bottom panel of Fig. 20. The Cp are the best-fit 
values for 155 < £ < 504, i.e., the weighted average of Ci 
in that multipole range, reported in Table II, and each 
data point is located at the logarithmic center of the en- 
ergy bin. The results of fitting Cp/(AE ) 2 oc E~ Tl are 
given in Table IV. Identifying Fi = 2r s , the best fit of the 
energy dependence suggests that the anisotropy is con- 
tributed by a source class with a power-law photon spec- 
trum characterized by F s = 2.4040.07 (2.33±0.08 for the 
cleaned data), assuming only one source class contributes 
appreciably to the anisotropy. As the single power-law 
energy dependence provides a very good fit to the data, 
attributing the anisotropv to a single source class is a 
plausible interpretation. 


We note that the spectral index implied for the dom- 
inant source class contributing to the anisotropy is in 
excellent agreement with the mean intrinsic spectral 
index of blazars as inferred from the Fermi-detccted 
members [10], strongly supporting the interpretation of 
the measured anisotropy as originating from unresolved 
blazars. We caution, however, that due to the varia- 
tion between individual blazars’ spectral indices, as well 
as possible effects of EBL attenuation and redshifting, 
the fluctuation angular power from blazars could exhibit 
some energy dependence in the range considered here. 
Therefore, assuming that blazars are the dominant source 
class contributing the anisotropy could lead to tension 
with the flatness of the measured fluctuation anisotropy 
energy spectrum. Additional support for a blazar inter- 
pretation could be provided by a detailed study of the 
energy-dependent anisotropy arising from specific blazar 
population models, calibrated to match the properties of 
Fermi- detected blazars, and the consistency of the pre- 
dicted anisotropy of these models with the measured am- 
plitude of the angular power. We defer a careful treat- 
ment cf this subject to future work. 


VIII. DISCUSSION 


Prior work has generated predictions for the angular 
power spectra of several source populations which may 
contribute to the IGRB. In most cases, the predictions 
for the anisotropy of the emission from a single source 
class have been cast in terms of fluctuation angular power 
C//(/) 2 , where Ci is the intensity angular power spec- 
trum of the source class and (I) its mean collective inten- 
sity in a specified energy range. Since the intensity con- 
tributions of most gamma-ray source classes to the IGRB 
are subject to large uncertainties, it is convenient to con- 
sider the fluctuation angular power, since this quantity is 
independent of the overall normalization of the intensity. 
This convention is particularly useful when the spatial 
distribution, number density, and relative flux distribu- 
tion of the sources is known or modeled, and the uncer- 
tainty in the collective intensify can be translated into a 
multiplicative factor that uniformly scales the observed 
intensity in all sky directions. For this reason, the fluctu- 
ation angular power is very well suited for characterizing 
an indirect dark matter signal since the intensity nor- 
malization scales linearly with the assumed annihilation 
cross-section or decay rate. 

By comparing the measured fluctuation angular power 
with predictions for various source classes, we can place 
constraints on the fractional contribution from each 
source class to the total intensity by requiring that the 
fluctuation angular power of the total emission is not ex- 
ceeded. Assuming that each contributing source class is 
uncorrelated with the others and the Poisson component 
dominates the angular power spectrum of each source 
class at the multipoles considered, the intensity angular 
pov/er of the total emission is given by 


Cp,tot = Cp i 4- Cp } 2 4 • •• (12) 


and so the fluctuation angular power of the total intensity 
is 

= c p,i . c p. 2 , fl3 \ 

{ho,? ( /to t) 2 (Act) 2 •" { > 

Rewriting the fractional contribution from source class i 
to the total intensity /, = (7i)/(/ to t), 


Cp ? tot 

(/tot) 2 “ h (/l) 2 /2 </ 2 ) 2 


(14) 


If we allow a single source class i to contribute all of the 
measured angular power, the source class is constrained 
such that 


f 2 ^ gp,tot/(/tot) 2 

Ji ~ Cp,i/(h) 2 ' 


(15) 


Source classes whose predicted fluctuation angular power 
exceeds the measured fluctuation angular power therefore 
cannot contribute the entirety of the measured intensity. 

It is important to note that the total intensity of the 
IGRB in this analysis is not equivalent to the intensity of 
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the isotropic emission reported in [5] , since that analysis 
employed much more stringent selection cuts to remove 
charged particle contamination, and used a fitting pro- 
cedure to remove contributions from non-isotropic com- 
ponents. In this analysis the total intensity of the emis- 
sion is simply the intensity that remains after the mask 
is applied, which may include some emission from non- 
isotropic components, as well as a non- negligible amount 
of charged particle contamination. However, we empha- 
size that since the charged particle contamination is pre- 
sumed to be nearly isotropic, with any potential fluc- 
tuations confined to large angular scales, it should not 
contribute to the intensity angular power at the multi- 
poles considered here (l > 155), and so a more robust 
comparison of models with the data could be achieved 
by comparing the predicted intensity angular power to 
the measurement. 

We now compare our measurement to existing predic- 
tions from the literature for the angular power spectra of 
various gamma-ray source classes, and summarize these 
results in Table V. We caution that the predicted angu- 
lar power can depend sensitively upon the adopted source 
model (in particular the shape of the flux distribution) , 
the assumed source detection threshold, and, for cosmo- 
logical source classes, assumptions regarding the effect 
on the observed energy spectrum of attenuation of high- 
energy’ photons by interactions with the EBL. Conse- 
quently, the constraints derived in this section should be 
taken only as indicative values for these source popula- 
tions. 

Ref. [27] predicted the fluctuation anisotropy from un- 
resolved blazars Cp/(/) 2 ~ 2 x 10“ 4 sr at £ ~ 100 (see 
Fig. 4 of that work). This is a factor of ~ 20 larger than 
the fluctuation angular power of ~ 10” 5 sr measured in 
the data, which suggests that emission from blazars, as- 
suming the model adopted in that study, contributes less 
than ~ i/y/20 ~ 20% of the total intensity. Note, how- 
ever, tfiat the ffux threshold for sources in the 1FGL 
catalog is between 0.5 and 1 x 10“ 9 photons cm' 2 s' 1 
for |6| > 30°, higher than the threshold assumed in [27]. 
If the blazar luminosity function is identical to the one 
assumed in [27], this discrepancy in thresholds would im- 
ply that the prediction for the blazar anisotropy in [27] is 
underestimated with respect to the one applicable to our 
analysis, since our masked maps include more bright un- 
resolved blazars. As a result, the constraint on the frac- 
tional intensity contribution to the IGRB from blazars 
for this model from our measurement would, if anything, 
be stronger. 

In contrast to the larger anisotropy expected from 
blazars, the fluctuation angular power at t ~ 100 pre- 
dicted for star-forming galaxies by Ref. [30] is ~ 2 x 
10 -7 sr at 1 GeV, far below the value measured in this 
analysis. Since star- forming galaxies would thus pro- 
vide a subdominant contribution to the measured angular 
power, this anisotropy measurement does not constrain 
their contribution to the total IGRB intensity. 

The anisotropy from dark matter annihilation in ex- 


tragalactic structures is predicted to be slightly smaller 
than that from unresolved blazars, although estimates 
can vary substantially due to differences in the adopted 
models. Moreover, for extragalactic dark matter anni- 
hilation the amplitude of the expected anisotropy can 
be highly sensitive to the energy spectrum of the emis- 
sion. The source energy spectrum depends on the dark 
matter particle mass and dominant annihilation chan- 
nels, while the observed energy spectrum is affected by 
redshifting and EBL attenuation. These factors can in- 
troduce a non-trivial energy dependence into the am- 
plitude of the anisotropy, particularly for high mass 
(~ 1 TeV) dark matter candidates. As a benchmark 
range, Refs. [23, 27, 39] predict the anisotropy from an- 
nihilation of extragalactic dark matter to be ~ 10~ 6 - 
10 -5 sr at i ^ 100 at energies of a few GeV, comparable 
to the measured value. 

The anisotropy from annihilation in Galactic dark mat- 
ter substructure is expected to be much larger than that 
from extragalactic dark matter. While variations in the 
assumed properties of Galactic substructure can lead to 
order-of-magnitude or larger variations in the predicted 
angular power, for typical assumptions the predicted fluc- 
tuation angular power is ~ 5 x 10” 5 sr at t *** 100 (e.g., 
Model A1 in Ref. [33]), which implies that dark matter 
annihilation can contribute less than ~ 43% of the total 
intensity. However, adopting alternative models for the 
substructure properties can increase or decrease the pre- 
dicted angular power by as much as ~ 2 orders of magni- 
tude [32-34], so the measured angular power represents 
a strong constraint on some substructure models. 

Galactic gamma-ray millisecond pulsars (MSPs) have 
also been considered as possible contributors to the in- 
tensity and anisotropy of the IGRB due to their extended 
latitude distribution [15, 31]. The emission from Galac- 
tic MSPs is expected to feature very large fluctuation 
anisotropy due to the relatively low number density of 
this source class compared to dark matter substructure 
or extragalactic source populations. Ref. [31] predicts 
fluctuation angular power at high Galactic latitudes of 
~ 0.03 sr at £ ~ 100 for this Galactic source class, which 
implies a contribution to the total IGRB intensity of no 
more than a few percent. 

In addition to the specific source populations consid- 
ered in this section, other Galactic source populations for 
which anisotropy predictions do not yet exist in the lit- 
erature may also contribute to the anisotropy as well as 
the intensity of the high-latitude diffuse emission. These 
include normal pulsars, as well as populations currently 
too faint to have had individual members detected by 
Fermi. The properties of these populations can be con- 
strained by both low-latitude and high-latitude source 
count analysis (in the case that individual members have 
been detected) [55], and also by the anisotropy analysis 
described in this study. We leave the detailed study of 
this to future work. 

We note that constraints derived in this section have 
not taken into account information about the likely en- 
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TABLE V: Maximum fractional contribution of various source populations to the IGRB intensity that is compatible with 
the best- fit constant value of the measured fluctuation angular power in all energy bins, (Cp/(J) 2 ) = 9.05 X 10 6 sr for the 
default data analysis or ( Cp/(7 ) 2 ) = 6.94 x 10“ 6 sr for the Galactic-foreground-cleaned data analysis. Indicative values for the 
fluctuation angular power Ct/(I) 2 of each source class are taken from existing literature (see text for details) and evaluated at 


i = 103. 


Source class 

Predicted Cioo/(I ) 2 

Maximum fraction of IGRB intensity 


M 

DATA 

DATA : CLEANED 

Blazars 

2 x 1(T 4 

21% 

19% 

Star-forming galaxies 

2 x HT 7 

100% 

100% 

Extragalactic dark matter annihilation 

1 x 1(T 8 

95% 

83% 

Galactic dark matter annihilation 

5 x 1(T 8 

43% 

37% 

Millis 3 Cond pulsars 

3 x 1CT 2 

1.7% 

1.5% 


ergy spectrum of the dominant contributing population, 
discussed in §V11, which is incompatible with sources 
known or expected to feature spectral peaks at the ener- 
gies we consider (for example, Galactic and extragalac- 
tic dark matter and MSPs). A careful study combining 
all observables obtained in this work would almost cer- 
tainly yield stronger constraints on contributing popula- 
tions. Furthermore, we have discussed the constraints 
obtainable on specific source populations by requiring 
that the total anisotropy from each population does not 
exceed the measured value. We emphasize, however, 
that stronger bounds could be derived if some fraction 
of the total anisotropy could be robustly attributed to 
one or more confirmed source classes, thereby reducing 
the anisotropy available to additional contributors. 


IX. CONCLUSIONS 

The statistical properties of the IGRB encode detailed 
information about the origin of this emission. The ad- 
vanced capabilities of the Fermi LAT. most notably its 
improved angular resolution and large effective area, have 
enabled a sensitive measurement of small angular-scale 
anisotropies in the IGRB. Using ~ 22 months of data, 
we performed an angular power spectrum analysis of 
the high-latitude diffuse emission measured by the Fermi 
LAT. Significant angular power above the photon noise 
level is detected in the data at multipoles 155 < i < 504 
in three energy bins spanning 1-10 GeV, and is measured 
at lower significance in the 10-50 GeV energy bin. The 
primary limitation of the measurement at high energies 
is low event statistics, which results in the measurement 
uncertainties being dominated by the photon noise. In 
this regime the measurement uncertainties scale roughly 
inversely to the number of events, and hence increasing 
the statistics by a factor of 2 or 3 could lead to a large 
enough improvement in the sensitivity of the analysis to 
allow a confident detection of angular power in this en- 
ergy range and greater sensitivity to energy-dependent 
anisotropy. 

The angular power measured in the data at 155 < £ < 


504 is consistent with a constant value within each energy 
bin, and the scale independence of the signal suggests 
that it originates from one or more unclustered popula- 
tions of point sources. Comparing the measured angular 
power with predictions for known and proposed gamma- 
ray source classes, constraints can be obtained on the 
collective intensity and properties of source populations 
that contribute to the IGRB. The fluctuation angular 
power detected in this analysis falls below the level pre- 
dicted for many source classes, including blazars, MSPs, 
and some scenarios for dark matter annihilation in Galac- 
tic and extragalactic structures. In these cases the mea- 
sured amplitude of the fluctuation angular power lim- 
its the contribution to the total IGRB intensity of each 
source class. 

The measured fluctuation angular power is consistent 
with a constant value over the energy range considered, 
however, due to the relatively large measurement uncer- 
tainties and limited number of energy bins, a mild energy 
dependence in this quantity cannot be excluded. The 
absence of a strong energy dependence in the fluctua- 
tion anisotropy energy spectrum suggests that a single 
source class may provide the dominant contribution to 
the anisotropy while providing a constant fractional con- 
tribution to the intensity of the IGRB over the energy 
range considered. We caution, however, that this anal- 
ysis is not sensitive to structure in the anisotropy en- 
ergy spectrum that is confined to small energy ranges, 
since the requirement of large event statistics to detect 
anisotropies at the measured level precludes fine energy 
binning of the data. We anticipate that future analyses 
that draw on larger data sets will be more sensitive to 
localized features in the anisotropy energy spectrum. 

The energy dependence of the intensity angular power 
of the data is well-described by that arising from a single 
source class with a power-law photon spectrum with in- 
dex T s = 2.40 ± 0.07. Interestingly, this value closely 
matches the mean intrinsic spectral index for blazars 
as determined from recent Fermi LAT measurements. 
While alternative scenarios invoking contributions from 
more than one source class to explain the energy depen- 
dence of the angular power are in principle possible, the 


27 


interpretation of the measured power as originating from 
a single source class with a power-law energy spectrum 
is an excellent fit to the data. To identify a specific 
population or populations as the source of the measured 
IGRB anisotropy, detailed analysis of population models 
for plausible source classes will be essential in order to 
verify that both the predicted intensity energy spectrum 
of the IGRB and the corresponding anisotropy signal pro- 
vides a consistent explanation of the data. 
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